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Abstract 

This paper proposes a probabilistic approach for acoustic emission (AE) source localization in 
isotropic plate-like structures based on an extended Kalman filter (EKF). The proposed 
approach consists of two main stages. During the first stage, time-of-flight (TOF) 
measurements of Lamb waves are carried out by a continuous wavelet transform (CWT), 
accounting for systematic errors due to the Heisenberg uncertainty; the second stage uses an 
EKF to iteratively estimate the AE source location and the wave velocity. The advantages of 
the proposed algorithm over the traditional methods include the capability of: (1) taking into 
account uncertainties in TOF measurements and wave velocity and (2) efficiently fusing 
multi-sensor data to perform AE source localization. The performance of the proposed 
approach is validated through pencil-lead breaks performed on an aluminum plate at 
systematic grid locations. The plate was instrumented with an array of four piezoelectric 
transducers in two different configurations. 

(Some figures may appear in colour only in the online journal) 



1. Introduction (i.e. extensional modes), whilst antisymmetric modes have 

displacements transverse to the wave propagation direction 

In the past two decades, significant effort has been made (i.e. flexural modes). Damage location using Lamb waves, 

towards the development of structural health monitoring can be achieved by either an 'active-passive' approach or by 

(SHM) systems based on guided ultrasonic waves (GUWs) a 'passive-only' approach. In the 'active-passive' approach, 

in order to increase safety and performance of aircraft diffractions of piezoelectrically actuated waves can be used to 

systems [1-3]. Damage events, both on ground and in-flight, locate existing damage in a post-mortem mode [4, 5]. In the 

can represent a threat to aircraft systems. Sources of 'passive-only' approach, growing damage or sudden impacts 

damage include hail impact, lightning strike, transport can be located by monitoring acoustic emissions in a real-time 

and handling, and foreign objects. Techniques based on mode [6-12]. This paper, focus on a 'passive-only' approach 

sparse arrays of sensors, which have the capability of using a sparse array of piezoelectric transducers, 

transmitting and receiving Lamb waves are among the Traditionally, damage or impact location is based on 

most promising candidates. The application of Lamb waves time-of-flight (TOF) triangulation of wave measurements 

can be challenging as they are multimode (many vibrating taken at multiple receiving points. This technique works very 

modes can propagate simultaneously) and dispersive (the well when the wave velocity (V g ) in the test material and 

propagation velocity depends on the wavefrequency /). In the arrival time (t{) of the signal at all three sensor locations 

plate-like structures two fundamental modes can propagate: are known. The damage or impact location is identified by 

symmetric (S) and antisymmetric (A). Displacements of the drawing three circles of radii (/?/), whose centers coincide 

symmetric modes occur in the direction of wave propagation with the three sensor locations. The radius R( is obtained 
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Figure 2. AE source location and sparse array of fz piezoelectric 
sensors. 
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Figure 1. AE source localization using triangulation (a) without 
uncertainty and (b) with uncertainty. 



by multiplying the time of arrival of the signal (t{) with 
the wave velocity (V g ). The intersection point of these three 
circles is the damage location (figure 1(a)). However, TOF 
and wave velocity are two uncertain parameters. In general, 
uncertainties can be caused by random and systematic errors. 
The random errors are caused by unknown and unpredictable 
changes in the TOF measurements, including instrumentation 
noise, temperature changes etc [13-18]. Systematic errors 
are mostly caused by the digital signal processing technique 
used for analyzing the time waveforms [19]. As a result, 
rather than the damage being located at a single point at 
the intersection of the circles as in figure 1(a), the damage 
can be located anywhere in the dark overlapped region in 
figure 1(b). Model-based approaches [20] have been proposed 
for impact location and identification in simple structures. 
In these methods the location and type of impact in a 
model are iteratively changed until the predicted responses 
match the measured responses. For damage/impact location 
in more complex structures, where accurate model-based 
predictions may not be possible, artificial neural networks 
have been proposed [21, 22]. However, neural networks 
require an extensive number of training observations prior 
to deployment, making these methods quite onerous from 
a computational or data storage point of view. Kundu et al 
[7], proposed an alternative method based on optimizing an 
objective function to find the impact location in isotropic and 
anisotropic plate. However, the proposed objective function 
in that reference had the inherent problem of multiple 



singularities which was overcome in Kundu et al [23] by 
modifying the objective function. The optimization technique 
was further improved by Hajzargerbashi et al [24]. In [25] 
and [26] the nonlinear least square optimization adopting 
the Gauss-Newton method was proposed to determine the 
location, time lag, and velocity of 'synthetic' AE signals. 
Ciampa and Meo [12] used Newton's iterative method to 
calculate the coordinates of the impact location and the 
wave velocity. A genetic algorithm (GA) for optimization 
of an objective function based on the modified triangulation 
methodology was proposed in [27]. 

Here an alternative probabilistic approach, based on 
the extended Kalman filter (EKF) theory, is proposed for 
AE source localization and wave velocity identification 
in isotropic plate-like structures. This approach, in which 
TOF and wave velocity are considered as Gaussian random 
variables, consists of two main stages. During the first 
stage, TOF measurements of Lamb waves are carried out 
by a continuous wavelet transform (CWT) accounting for 
systematic errors due to the Heisenberg uncertainty; the 
second stage uses the EKF to iteratively estimate the location 
of the AE source and the wave velocity. The paper begins 
by briefly introducing the algorithm for the AE source 
localization. Section 3 summarizes the main characteristics 
of the CWT for TOF and frequency measurements. This 
is followed by section 4 which briefly introduces the EKF 
theory. In section 5, the AE source location based on the EKF 
shall be discussed. Section 6 reports the experimental setup 
followed by the results in section 7. Finally, the conclusions 
are given in section 8. 

2. Source location algorithm 

Assuming an arbitrary Cartesian coordinate system, the AE 
source is at unknown coordinates (x s , v s ) in the plane of 
the plate and the sensors are located at (x/, v/), as shown 
in figure 2. It is well known that given the wave velocity 
of a particular Lamb wave mode, at least three sensors are 
necessary to locate an AE source in plate-like structures. If t m 
is the travel time to reach the first sensor (master sensor), and 
Atmi are the time difference between master sensor and the /th 
sensor, the following equations can be obtained: 



(x m - x s y + (y m - y s y - (t m v e y = o 



(2.1) 
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(xi - x s ) 2 + (y t - y s ) 2 - [(t m + At mi )V g ] 2 = (2.2) 

where V g is the group velocity of the particular Lamb wave 
mode, which is a function of the product of the frequency / 
and the plate thickness d: 



V g = Fifd). 



(2.3) 



Knowing the properties of the plate, mass density p, thickness 
d, Young's modulus E, and Poisson ratio v, V g can be 
calculated by the Rayleigh-Lamb equations [28]. Combining 
equations (2.1) and (2.2) the following equation can be 
derived: 



At n 



(x/-x s ) 2 + (v/-vs) 2 - 



(x m -X s ) 2 + (Vm- V S ) 2 



-. (2.4) 



In general, the differences between the TOF of the master 
sensor to different sensors (Af m ;) are used to solve this 
set of nonlinear equations with unknowns X = (x s , v s , Vg). 
However, as mentioned in section 1 , several sources of error 
in the TOF measurements (i.e. dispersion, noise, temperature, 
etc), can affect the accuracy of this solution. In this work, to 
take into account uncertainties in the arrival time or TOF t[ 
and the wave velocity V g , the arrival time t\ and unknowns x s , 
v s , V g are treated as mutually independent Gaussian random 
variables. Based on this assumption, the probability density 
function of the time difference Af m ; can also be defined as 
a Gaussian random variable [29] with mean and variance 
defined as: 



At n 



r t+°t 



(2.5) 



Time differences At m t are determined through a 
time-frequency analysis using a CWT which has inherent 
uncertainties in time and frequency localization caused by the 
Heisenberg inequality, as will be explained in section 3. In 
this probabilistic framework, an extended Kalman filter (EKF) 
approach is proposed in section 4 for the optimal estimation 
of the state vector X = (x s , v s , V g ). 

3. The continuous wavelet transform (CWT) 

A variety of methods for increasing the accuracy of TOF 
measurements of acoustic waves have been developed, includ- 
ing the short time Fourier transform (STFT), Wigner-Ville 
distribution (WVD), CWT, and Hilbert-Huang transform 
(HHT) [3]. In this work a CWT has been employed for its 
good resolution both in the time and frequency domain and 
for its multi-resolution property [30]. The CWT of a signal 
fit) is defined as [31]: 



WT(,s, b) 



1 ( + °° 

V s J -co 



fm* 



dt 



(3.1) 



where s > is the scale or dilatation parameter, b is the 
translation parameter localizing the wavelet in the time 
domain. ty(t) is called the mother wavelet, and ty*(f) 
denotes the complex conjugate of the mother wavelet. Several 
researchers have used the CWT to analyze Lamb waves in 
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Figure 3. Time-frequency localization. 



plate-like structures for AE source localization [11, 12, 25, 
32]. However, in this paper a new approach is proposed to 
take into account systematic uncertainties, in the arrival time 
measurements and group wave velocity assumption, due to the 
Heisenberg uncertainty principle. 

3.1. Heisenberg uncertainty 

The local resolution of the CWT in time and frequency 
domains is limited according to the Heisenberg inequality [31] 
which is defined as: 



^22 ^ 1 



(3.2) 



where 07 and a M , representing the local time and frequency 
resolution of the CWT, depend on the scale parameter is) and 
on the mother wavelet properties as: 



a t = o t ^s, 



&coV 



(3.3) 



Here 07^ and o^ are the standard deviations representing 
time and frequency resolution of the mother wavelet, defined 
as [31]: 



1 



&M 



l^ll 



/ 



+00 



(r-r*) 2 |*(0l 2 dr (3.4) 



CTffl* 



11*112 V -/— 00 



icD-coy) 2 \i>icD)\ 2 dcD (3.5) 



where % and coy are the centers of ^it) and ^ico) 
respectively and can be calculated by: 



ty 



coy 



I 



+°° JV(Q| 2 

-00 11*11! 



/ 



+°° \V(CD)\ 2 A 

co y~ dco. 

Il^llo 



(3.6) 



(3.7) 



Here || • H2, and | • | denote the classical norm two and 
the absolute value respectively. According to equation (3.2) 
an improvement in time localization causes the deterioration 
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in frequency localization and vice versa. For some mother 
wavelets, such as the complex Morlet, the localization in 
the time and frequency domain can be represented by two 
Gaussian probability density functions with variance <r r 2 
and <t 2 and mean values of b + sty and ^f-, respectively. 
Figure 3 shows these probability density functions and their 
corresponding Heisenberg boxes of two wavelets ^bi,si(0 
and ^b 2 ,s 2 (0- In this work, the variance in time domain <r r 2 
will be used to define the uncertainty in arrival time, and 
the variance in the frequency domain a 2 will be correlated 
to the uncertainty in the group wave velocity defined by the 
variance <jy . It can be stated that by increasing the time 
resolution the systematic error in arrival time determination 
decreases. On the other hand, the uncertainty in the initiation 
of group velocity <jy , which is related to a 2 through the 
well-known Rayleigh-Lamb equations, increases according 
to equation (3.2). 

3.2. Complex Morlet wavelet 

There are several mother wavelets reported in the literature 
and exploited to extract time and frequency information 
from acquired signals. In this study a complex Morlet was 
employed as it is able to separate amplitude and phase 
information, but also it achieves a desirable compromise 
between time and frequency resolution, with the area of the 
Heisenberg box equal to 0.5 [31]. The complex Morlet is 
defined as: 



*(0 



1 



=e ^e" 



2nf c jt 



(3.8) 



where f c and/b represents the center frequency and frequency 
bandwidth, respectively, fi, controls the shape of the mother 
wavelet and, for a given scale parameter s, it determines the 
resolution in the time and frequency domain. The following 
equations show the a t and a M in equation (3.3) for the complex 
Morlet [33]: 



Of = 



sy/jb 



1 



sVJb 



(3.9) 



Figure 4 shows the real and imaginary part of the complex 
Morlet wavelet and the corresponding Gaussian distribution 
functions with mean t\&= and standard deviation a t y — ^ 
in the time domain and with mean coy = 27tf c and standard 
deviation o^ — -j= in the frequency domain. To find the 
arrival times a scale s should be selected, as explained in 
section 3.3. 
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Figure 4. Complex Morlet wavelet (a) real part in the time domain, 
(b) imaginary part in time domain, (c) probability distribution 
function in the time domain and (d) probability distribution function 
in the frequency domain. 



3.3. Determining the difference of arrival time 

The squared modulus of the CWT (scalogram) represents the 
energy density of a signal in the frequency-time domain [31]. 
The scalogram can be defined as: 



|WT(5, b)\ 2 = WT(s, b) • WT*(s, b) 



(3.10) 



where WT* represents the conjugate of wavelet transform 
WT. The maximum value of the coefficients of the scalogram 



is achieved at the instantaneous frequency [11, 25, 30], 
corresponding to the dominant frequency in the signal 
analyzed at each instant in time. These coefficients, taken 
at the instantaneous frequency in a time-frequency domain, 
determine the ridges [31]. The projection of the ridge 
corresponding to the dominant frequency ft represents the 
arrival time of the wavepacket. In this work, because the 
dominant frequency can be different at each sensor, the 
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average/ is used to calculate the arrival time tf. 

- f= m. 



(3.11) 



(3.12) 



The corresponding scale s is defined as: 

&)* 

where coy is defined in equation (3.7) and is co^ — 2nf c for 
the complex Morlet wavelet. By substituting equation (3.12) 
into (3.9), cr t and a M for the complex Morlet wavelet can be 
defined as: 



crt = 



2/ 
/ 

fc+Jfb 



(3.13) 



(3.14) 



Substituting equation (3.13) into (2.5) the variance of the 



difference in arrival times a\ f is defined as: 



2 _ r > ( U4K \ 



(3.15) 



This value is used to model the uncertainty of the 
measurement vector, which is given as input to the proposed 
EKF algorithm. For the sake of clarity a brief general 
derivation of the EKF theory is presented in section 4. 

4. Extended Kalman filter (EKF) 

In general, a Kalman filter (KF) is a recursive data processing 
algorithm that estimates the state of a noisy dynamic 
system [34]. When we talk about the state of a system we 



mean a vector X consisting of n variables describing some 
interesting properties of the system. To estimate the state 
a KF processes all available measurements, both accurate 
and inaccurate measurements. The KF has been widely 
used in a wide range of applications due to its simplicity, 
optimality, robustness, sensor fusion ability and its capability 
in considering the uncertainties due to the mathematical 
model and measurements [34-37]. In general the KF is 
a two steps process: (1) state prediction according to a 
mathematical model and (2) correction of the state according 
to the measurements collected by sensors. More specifically, 
in the prediction step the KF estimates the state of the system 
at a given time instant and then obtains a feedback control 
in the correction step, by incorporating a new measurement 
result into the a priori estimate, in order to gain an improved 
a posteriori estimate. Although the underlying approach is 
very promising there is a critical limitation. In fact the 
KF assumes that system and measurements are adequately 
modeled by a linear dynamic system. The EKF is similar 
to the KF but it can be used in nonlinear systems. For this 
reason the EKF has been used to handle the set of nonlinear 
equation (2.4). 

4.1. EKF formulations 

Without loss of generality, let us consider a generic nonlinear 
dynamic system with no external inputs, as follows: 



X M =f(X k )+w k 
Z k = h(X k ) + v k 



(4.1) 
(4.2) 



where X^ is the current ^-dimensional state vector of 
the system, Z^ is the m-dimensional vector of current 
measurements at the k\h step, / and h are known functions, 
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and Vk ~ A/"(0, /?) and u^ ~ A/"(0, g) are uncorrelated noisy 
Gaussian processes with zero mean and covariance matrices R 
and Q. The objective is to obtain estimates Xk of the state Xk 

A T 

using measurements Z# so that tr(E[(Z^ — Z^)(Z^ — X&) ]) 
is minimized (tr, E, and T denote trace, expected value, and 
transpose respectively). According to the EKF, the / and h 
functions can be approximated through the first-order Taylor 
series expansions with respect to X. With this assumption, the 
prediction stage of the process may be analyzed through the 
following equations: 

(4.3) 
(4.4) 
(4.5) 



(4.6) 



where <$>k is the state transition matrix. X^, P^, and Z^ 
are the a priori (predicted) estimates of the state vector, 
the error covariance matrix, and the predicted measurements, 
respectively. The correction stage equations are expressed as: 



x k =f(X 


f-l) 


*t =**-i*t-i*it-i+Gt 


z k = h(X k ) 


% = df(x) 

dX 


x =K-i 



X+=X k +K k {Z k -Z k ) 

p+ = (/ - k k H k )p- 



(4.7) 
(4.8) 



where / is an identity matrix; X^ and P^, are a posterior 
(corrected) estimates and the Jacobian matrix H k contains the 



partial derivatives of the measurement function h with respect 
to the state X, evaluated at the prior state estimate X k and is 
defined as: 



H k 



dh{X) 



dx 



X=X7 



and Kk is the Kalman gain defined as follows: 

-l 



K k -- 



p-H T k (H k p-H T k +R) 



(4.9) 



(4.10) 



The EKF has been shown to be successful in many practical 
nonlinear applications [34]. To the best of the authors' 
knowledge, none of the techniques presented in the literature 
for AE source location take advantage of the EKF theory. 

5. AE source localization based on EKF 

In this work, the state of the system X consists of 
three parameters, the AE source location (x s , v s ) and the 
wave velocity V g . Given a sparse array of n sensors, 
the time differences Ar m ; between the /th sensor and the 
master sensor represent the measurement data vector Z = 
[Afm/]( n _i) X i, n > 3. The state of the system X is related to 
the measurement vector Z through the nonlinear measurement 
function h defined in equation (2.4). The key idea underlying 
the proposed approach is that the AE source location can 
be considered as an unmovable (static) or slowly fluctuating 
point. Under this assumption, the prediction step of the EKF 
can be simplified, as the transition matrix <&k reduces itself 
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Figure 7. Proposed approach flowchart. 




Figure 8. Experimental setup. 



to the unity matrix /. Furthermore, since the state X has 
no superimposed noise, Q is considered to be null. Under 
these assumptions, the EKF algorithm can be simplified 
by neglecting the prediction stage. Thus the posterior 
estimation is equal to the a priori estimation at any step k, 
(i.e. the superscripts '— ' and '+' can be neglected) and the 
equations (4.5), (4.7), (4.8) and (4.9) take the following forms: 



Z k = h(X k -i) 
X k = X k - 1 +K k (Z-Z k ) 

P k = (i-k k H k )P k - l 



(5.1) 

(5.2) 
(5.3) 



H k 



dh{X) 



dx 



X=Xk-i 



Also equation (4.10) can be rewritten as: 



K k : 



:P k -iHj(H k P k -iHj+Ry 



(5.4) 



(5.5) 



It should be noted that the measurement vector Z in 
equation (5.2) does not have the subscript k because its (n — 1) 
components (difference in arrival times) can be considered 
constant, once the AE signals have been acquired. Thus during 
the iteration to estimate the state vector X this measurement 
vector is constant. The uncertainty of measurement vector Z is 
modeled as a zero mean white Gaussian noise with covariance 
matrix R with diagonal terms defined in equation (3.15). 



R = 











^2 



m(«— 1) _ 



(5.6) 



The objective is to estimate the state vector X — [x s ,y s , V g ] 
given the measurement vector Z and covariance matrix R. The 
EKF will iteratively correct a priori knowledge of the state 
vector X and covariance matrix P, with respect to the Z and R. 
The iteration begins with the initialization of the state 
vector estimate Xq and its covariance matrix Pq: 



*o = [*so,5v^g ] 



/\) = 



- XsO 



o o 



<o 
< 







(5.7) 
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Figure 9. Initial Gaussian locations in the first configuration of sensors and actual locations of simulated AE sources: (a) 3D view and 
(b) planar view. 



5.7. Initiation of location 

The starting estimates of the AE source location x So and % 
can be calculated by geometrical considerations. Given an 
arbitrary group of ft-fired sensors at locations (xi,yf), it can 
be assumed that estimated coordinates of the source location 
x s and y s can take any value in the intervals of (x s l, x s tj) and 
(y S L, .Vsu)? where x s l and j s l are the lower bounds defined as: 



x sL = min(Mn); JsL = min(£yj]„) 

and x s u and y s \j are the upper bounds defined as: 



(5.8) 



x sU = max([*;]„); 



y s u = max (tyi]„). (5.9) 



For the sake of clarity in figure 5 is shown the simplest 
case of three fired sensors. Making these assumptions, mean 



and variance for an uncorrelated uniform distribution can be 
expressed as: 



(x s u + * s l) 



^s 



(VsU + JsL) 



(5.10) 



% = 12 ( X sU _ X *l) 



y*o 



12 



(jsU-JsL) • 



(5.11) 



5.2. Initiation of group velocity V g 

In section 3 the frequency if) has been assumed as a 
Gaussian random variable with mean and variance defined 
in equations (3.11) and (3.14). The starting estimates of 
the group wave velocity V g0 and its variance oVg can 
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Figure 10. Initial Gaussian locations used in the second configuration of sensors and actual locations of generated AE sources: (a) 3D view 
and (b) planar view. 



be calculated determining the relationship between the 
probability distributions of the two random variables/ and V g 
when they are related by V g = F(fd). Assuming the functional 
form F{fd) is given and deterministic, we are interested in 
determining the probability density function of V g given the 
probability density function off (the thickness of the plate d 
is deterministic). Because in general the functional form of 
F is not available, in this work an unscented transformation 
(UT) has been used. The UT is a method for calculating the 
statistics of a random variable which undergoes a nonlinear 
transformation [38]. Assuming the random variable fd has 
mean fd and variance oh, three weighted sigma points xo> 
XI, and xi are necessary to calculate the statistics of V g . These 



points are given by: 

Xo=fd, Wo 



Xl 



(1+K) 

W\ = 



X2 



w 2 = 



2(1 + k) 
1 

2(1 +k) 



(5.12) 



where Wt is the weight associated with the iXh point and 
k is an arbitrary number providing l + /c^0;/c = 0is 
chosen in this work. Given the set of sigma points calculated 
by equation (5.12), the transformation of the Gaussian 
probability density function of the random variable fd to the 
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Figure 11. Each figure depicts the time histories of the four signals acquired by the sensors, the contour plot of the scalogram of the CWT, 
and the line profile of the scalogram showing the procedure to find the time of arrival at/ = 137.00 kHz for AE point 4. 



approximated Gaussian distribution of V g is summarized in (2) The mean of random variable V g can be calculated from: 
the following steps. 

(1) Calculate the group velocity corresponding to each 

sigma point according to dispersion curvature or the ^ The variance of random variable V g can be determined by: 

Rayleigh-Lamb equations. 2 

^^{^-VgH^-Vg} 1 . (5.15) 



Vi = F(xt), 



i = 0,1,2. 



(5.13) 



-i 
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Figure 12. AE source localization results using three sensors in the 
first setup configuration. 



Table 1. Properties of four initial Gaussian locations. 



Initial Gaussian location I 



II 



III 



Note: The origin of Cartesian coordinates is at sensor 1. 

Table 2. Properties of initial velocity in the first setup 
configuration. 



IV 



x sU (cm) 


12.7 


25.4 


25.4 


-25.4 


Vsu (cm) 


12.7 


25.4 


-25.4 


25.4 


<*l (cm 2 ) 


13.44 


53.76 


53.76 


53.76 


tf s (cm 2 ) 


13.44 


53.76 


53.76 


53.76 


x SQ (cm) 


6.35 


12.7 


12.7 


-12.7 


Js (cm) 


6.35 


12.7 


-12.7 


12.7 



Actual location /(kHz) Vg (ms l ) 



V 9.i 



(ms" 1 ) 2 



1 


152.75 


2901.8 


1063.7 


2 


150.00 


2892.1 


1076.2 


3 


112.25 


2724.8 


1612.4 


4 


137.00 


2843.4 


1258.7 
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Figure 13. AE source localization results using four sensors in the 
first setup configuration. 



Now the initial group wave velocity for a specified 
mode is V g0 = V g and its corresponding variance used in 
equation (5.7) is <Xy . Figure 6 shows this transformation 

for a given frequency/ assuming Ao as the dominant Lamb 
wave mode. In the appendix is presented another approach 
that can be used to perform this transformation in the case 
where the functional form of F is known. Overall, the 
proposed approach consists of two stages as shown in the 
flowchart in figure 7; the first stage provides information 
about time of arrivals, frequency, and their uncertainties; the 
second stage uses a priori information about the states and 
iteratively estimates the state vector according to the measured 
information obtained from the first stage. 

6. Experimental setup 

Experiments were carried out using an aluminum plate with 
dimensions 910 mm x 910 mm x 3.175 mm to validate 
the proposed algorithm. The plate was instrumented with 
an array of four piezoelectric transducers in two square 



configurations, with dimensions 254 mm and 560 mm, 
respectively. The experimental setup is shown in figure 8. For 
the data acquisition, a LeCroy Waverunner oscilloscope with 
sampling frequency of 10 MHz was used. Each sensor was 
connected to a preamplifier. During testing, AE sources were 
generated by pencil-lead breaks at systematic grid locations. 
To preferentially generate the antisymmetric (Ao) Lamb wave 
mode the lead was fractured on the plate surface [39]. 
Post-processing of the received signals was performed with 
a PC running a Matlab software code implemented by the 
authors. The complex Morlet parameters were set by trial and 
error and fixed (f c = 3 Hz,/b = 0.7 Hz). 

The first configuration was used to investigate the 
convergence behavior of the EKF to the initial guess, Xo and 
Pn- In particular, four initial Gaussian locations were used, as 
shown in figure 9. In this figure the actual source location is 
marked as a squared dot (□) and each sensor with a circle 
(CO- The lower bound x s l and v s l is the location of the 
master sensor (sensor 1) located at the origin of the coordinate 
system and the upper bounds x s tj and v s u are summarized in 
table 1. The second configuration was used to validate the 
performance of the proposed approach in additional known 
locations. In particular, artificial AE sources were generated 
by pencil-lead breaks on 16 grid points. The initial Gaussian 
location was calculated based on equations (5.8) and (5.9) 
where the lower bound x s l and v s l is the location of the master 
sensor (sensor 1) located at the origin of the coordinate system 
and the upper bound x s tj and v s u is the location of sensor 3. 
Figure 10 shows the initial Gaussian location used to perform 
the source localization at each point. 

7. Results and discussion 

The frequency / of each simulated AE source, and the 
corresponding mean and variance of the initial group velocity 
computed using the UT, are summarized in table 2. For the 
sake of clarity, figure 1 1 shows the received signals, contour 
plot of the scalograms, and the line profile at the average 
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Figure 14. Source location error in the x and y coordinates using four sensors and considering as initial Gaussian locations: ((a), (b)) I; 
((c), (d)) II; ((e), (f)) III, and ((g), (h)) IV. 



frequency / = 137 kHz, for an AE source simulated at a 
given point (i.e. point #4) in the first sensor configuration 
setup. The transformation of the probability density function 
from the frequency domain to the group velocity domain has 
been shown in figure 6. The error of the estimated AE source 
location (x Sk , y Sk ) at the Mi step can be defined as: 



*** = ioo 



Sy sk = 100 



ksal 

Kl 



% 



% 



(7.1) 



(7.2) 



where (x sa , v Sa ) are the coordinates of the actual source 
location. Figure 12 shows the results of the proposed approach 
using data from an array of three sensors (i.e. 1, 2, 4). In 
particular, the actual source location is marked as a squared 
dot (□). The estimated location, using the proposed algorithm, 
is represented on the figure by the symbol (+). It can be 
observed that using three sensors may lead to inaccurate 
estimation of source locations. In fact, the average errors 
of four simulated AE sources after 10 iterations in x and y 



directions according to equations (7.1) and (7.2) are 36% and 
23%, respectively. Source location results obtained by adding 
the information provided by a fourth sensor, are illustrated in 
figure 13. In this case, the average errors after 10 iterations 
for x and y estimations are 4% and 5%, respectively. In 
general, the more information in hand the more accurate 
results can be obtained. Because of the matrix structure, the 
proposed method is very flexible in adding or removing the 
information from a network of sensors at the data processing 
stage. Furthermore, it is possible to consider different modes 
in the measurement vector. Indeed, if in addition to the 
antisymmetric Ao, the symmetric mode So is excited, the 
number of measurements increases from (n — 1) to 2(n — 1), 
with just one more unknown added to the state vector X, 
that is, the So group wave velocity. The convergence of the 
proposed approach using an array of four sensors is shown 
in figure 14. These results were obtained considering as a 
starting point the four initial Gaussian locations illustrated 
in figure 9. It can be observed that the proposed algorithm 
is capable of converging in just a few steps with an average 
error of less than 5% in both x and y coordinates, and in 
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Figure 14. (Continued.) 



Table 3. Actual and estimated source positions in the first setup 
configuration (with the first initial location). 

Point no: x s (mm) a x s (mm) y s (mm) a % (mm) V g (ms _1 ) 



1 


50.8 


47.01 


50.8 


49.7 


2948.1 


2 


101.6 


98.81 


50.8 


43.25 


3022.4 


3 


50.8 


48.61 


101.6 


100.5 


2755.7 


4 


101.6 


98.08 


101.6 


97.25 


2883 



Note: The origin of Cartesian coordinates is at sensor 1. 

the worst case the error is less than 7.5%. It should be noted 
from figure 14, that the first Gaussian initial location can be 
considered as the worst scenario. In fact, although the first 
initial Gaussian location has the smallest starting error (i.e. it 
is closer to the actual locations), since the variance is small, 
the contribution of the correction stage is reduced; as a result 
the initial state tends to resist any correction or update. The 
results using the first initial Gaussian location are summarized 
in table 3. In addition, figure 14 demonstrates the robustness 
of the proposed approach to converging to the optimal solution 
in just six steps, regardless of the initial starting location. 

Results obtained using the second configuration are 
shown in figure 15. The average errors in x and y are less 



than 5% and 7% respectively. In this setup sensors cover 
the area 4.8 times the first configuration and the maximum 
distance between sensors is more than 560 mm. Table 4 
summarizes, for each point, the results in terms of: actual 
location coordinates (x Sa , v Sa ), estimated location coordinates 
(x s , y s ), and estimated wave velocity Vg. 

Furthermore, an interesting feature of the proposed 
method is the capability to perform sensor and data fusion. 
Indeed, suppose that two techniques with different reliabilities 
are used during the first stage to perform time of arrival 
measurements, a decision should be made to combine the 
information according to their degree of reliability. The 
EKF is capable of fusing this information, so increasing the 
accuracy of the source location. 

8. Conclusion 

In this paper a probabilistic framework for locating acoustic 
emission sources is proposed based on the EKF theory. 
The proposed approach, in which TOF and wave velocity 
are considered as Gaussian random variables, consists of 
two main stages. During the first stage, TOF measurements 
of Lamb waves are carried out by a continuous wavelet 
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Figure 15. AE source localization results using four sensors in the 
second setup configuration. 



Table 4. Actual and estimated source positions in the second 
configuration setup. 



Point no: 


x s (mm) a 


x s (mm) 


J s (mm) a 


% (mm) 


^(ms" 1 ) 


1 


70 


50.53 


70 


51.25 


3130.4 


2 


210 


201.4 


70 


78.64 


3121.8 


3 


350 


358.7 


70 


79.16 


3094.8 


4 


490 


510.1 


70 


49.16 


3110.8 


5 


70 


76.73 


210 


199.6 


3124.5 


6 


210 


206.4 


210 


206.4 


3117.6 


7 


350 


347.1 


210 


208.5 


3102.4 


8 


490 


478.7 


210 


199.2 


3142 


9 


70 


79.83 


350 


358.5 


3129.2 


10 


210 


207.7 


350 


350.9 


3127.5 


11 


350 


351.5 


350 


351.5 


3124.2 


12 


490 


480.1 


350 


361.1 


3087.7 


13 


70 


44.05 


490 


516 


3089.7 


14 


210 


194.9 


490 


483.3 


3101.7 


15 


350 


359.7 


490 


480.8 


3147.1 


16 


410 


509.4 


490 


513.1 


3129.4 



Note: The origin of Cartesian coordinates is at sensor 1. 



transform (CWT) accounting for systematic errors due to 
the Heisenberg uncertainty; the second stage uses the EKF 
to iteratively estimate the source location and the wave 
velocity. The advantages of the proposed algorithm over 
the traditional methods include the capability of: (1) taking 
into account uncertainties in TOF measurements and wave 
velocity and (2) efficiently fusing multi-sensor data to perform 
an efficient AE source localization. The proposed algorithm 
is computationally efficient and can be used in real-time 
applications. It should be noted that the signal processing 
stage (CWT) can be replaced by other methods capable of 
providing probabilistic measurement data for the estimation 
stage. The performance of the proposed algorithm has 
been validated through pencil-lead breaks performed on an 
aluminum plate at systematic grid locations. The plate was 
instrumented with an array of four piezoelectric transducers. 



Appendix 

If the closed form solution of the group velocity as a 
function of frequency-thickness (equation (2.3)) is known, 
the probability density function of random variable V g can be 
stated as: 



-l, 



Pv g (V g )=Pf d [F-\V g )] 






(A.l) 



where P denotes the probability density function and d 
is the differential operator. Equation (A.l) is held if the 
function V g = F(fd) is continuous and strictly monotonic [29]. 
It is well known that the Rayleigh-Lamb equation can 
be approximated, for low frequency-thickness by the plate 
theory for the two fundamental Lamb wave modes (So and 
Ao) [40]. Under these conditions the phase velocities for So 
and Aq are expressed as: 



y ps = 



V **o = 



P (1 - vf 

E(2jrfd) 2 
12,0(1 - v 2 ) 



(A.2) 



(A.3) 



where />, E, and v are density, Young's modulus, and Poisson 
ratio respectively. The group velocity V g can be calculated 
from the phase velocity V p h according to the following 
equation [41]: 



Ve = V: 



P h 



V ph - (fd) 



dV T 



P h 



d(fd) 



(A.4) 



The approximate group wave velocity of So and Ao can be 
expressed as: 



V *s = 



P (1 " v) 2 



V gA =F(fd)=2 



EilitfdY 
\2p(\ - v 2 ) 



(A.5) 



(A.6) 



More accurate approximation of Ao group velocity as an 
analytic function of frequency-thickness can be found in [42]. 
From equation (A.5) it is evident that the approximated So is 
not a function of frequency-thickness so the transformation 
of the probability density function of random variable fd 
to the approximated group wave velocity of So is not 
meaningful [28]. On the other hand from equation (A.6) it 
is obvious that the approximated group velocity for Ao is a 
continuous, strictly monotonic, and increasing function of fd 
so equation (A. 1 ) is applicable in transferring the probability 
density function of random variable /d to V gA . Let us proceed 
with defining the inverse function of (A.6): 



F-HVgj*) =fd 



A J 



V 2 
8tt 



I2p(l - v 2 ) 



(A.7) 
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Figure A.l. Transformed probability density functions of random 
variable V gA using analytic, Monte Carlo, and unscented methods. 



The derivative of equation (A.4) is: 
V* 



^~\v gAo ) 



dV E 



A 



An 



lip {I- v 2 ) 



(A.! 



If the probability density function of the random variable 
fd (Pfdifd)) is Gaussian with mean fd and variance ok the 
Pfdifd) can be defined as: 



Pfd (fd) = 



1 



- ifd-fdf 



(A.9) 



(TfdWlTT 

After substituting equations (A.7) and (A.8) into (A.l) 
and considering the Pfdifd) stated in equation (A.9), the 
probability density function of V gA is stated as: 



p v 8Ao (V gAo ) = 



y gA 

An 



1 



J^m^ _ r \ 



°fd 



12p(l - v 2 ) 



(A. 10) 



The probability density function of the approximated Ao 
Lamb group wave velocity is derived in equation (A. 10). 
An alternative approach to perform this transformation is 
using a Monte Carlo simulation. However, although this is a 
reliable method, it is severely time consuming. To compare 
the performance of these three methods (analytic, Monte 
Carlo, and unscented transformation) data obtained from the 
AE source generated at the point 4 have been used, as 
shown in figure A.l. The Monte Carlo simulation has been 
carried out by randomly generating 200 000 samples from the 
random variable fd and the group wave velocity determined 
from equation (A.6). It can be observed that the unscented 



transformation method is able to capture the probability 
density function of random variable V gA just using three 
sigma points. 
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